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摘要 


f ESR, P 范 深 ! 


射电 天 文 图 像 去 卷 积 是 射电 天 文学 中 的 一 项 关键 数据 处 理 技术 , 其 主要 目标 是 去 除 天 空 图 像 中 由 天 文 观测 
仪器 引入 的 效应 , 从 而 复原 原始 的 天 空 图 像 . 射电 望远镜 阵列 采用 稀 芍 干涉 阵列 , 成 像 原理 与 光学 望远镜 有 


所 不 同 . 如 果 UV 空间 中 的 采样 点 不 足够 密集 , 将 会 导致 在 图 像 重 建 时 无 法 获得 足够 高 分 辩 率 的 信息 , 传统 


的 射电 天 文 图 像 重建 算法 未 能 根本 解决 UV 空间 欠 采 样 的 问题 . 本 文 基 于 压缩 感知 理论 框架 , 并 结合 射电 
种 新 的 射电 天 文 图 像 去 卷 积 算法 , 该 算法 将 脏 图 的 去 卷 积 过 程 转 化 
为 一 个 旨 在 求解 全 局 最 小 化 的 凸 优化 问题 , 即 基于 IUWT-CS 的 射电 干涉 图 像 重 建 算法 . 为 了 评估 该 算法 
究 使 用 了 射电 天 文学 仿真 软件 包 OSKAR 对 SKA1-low 射电 望远镜 阵列 进行 模拟 观测 , 并 


天 文 图 像 稀 玖 性 的 先 验 知识 , 研究 了 


的 重建 性 能 ， 


对 观测 得 到 的 点 源 和 扩展 射电 源 进行 了 去 卷 积 处 理 . 实验 结果 表明 , 与 HOGBOM-CLEAN, MS-CLEAN 


fu IUWT-FISTA 方法 相 比 , IUWT-CS 方法 显著 提高 了 射电 图 像 的 重建 质量 , 实现 了 更 加 精细 的 去 噪 和 复 


原 效 果 . 


关键 词 。 图 像 去 卷 积 , 压缩 感知 , 稀疏 表示 , AREY 


PACS: 95.75.Mn,07.05.Pj,95.75.Tv 


在 射电 望远镜 阵列 观测 的 过 程 中 , 由 于 数据 采 
样 的 不 均匀 和 不 完备 ， 原始 观测 数据 在 成 像 时 会 
产生 类 似 点 源 扩散 函数 (Point Spread Function, 
PSF) 的 旁 鸭 效应 , 同时 , 方向 独立 效应 (direction- 
independent effects, DIES)! 和 方向 依赖 效应 
(direction-dependent effects, DDES) P! 的 校准 也 


会 对 成 像 有 较 大 的 影响 , 为 了 改善 UV 覆盖 不 均 义 
以 及 校准 对 成 像 的 影响 , 通常 利用 测 得 的 可 见 度数 
据 和 PSF 进行 去 卷 积 处 理 . 在 现代 射电 干涉 测量 
技术 中 , 多 种 成 熟 的 去 卷 积 技术 已 用 于 重建 射电 图 
像 , 例如 广泛 使 用 的 HOGBOM-CLEAN 算法 站 和 
多 尺度 CLEAN(Multiscale Clean, MS-CLEAN) 算 
ik. CLEAN 算法 旨 在 从 模糊 且 受 噪声 干扰 的 观 
测 数 据 中 重建 清晰 点 源 的 位 置 和 强度 信息 , 其 处 理 


结果 通常 由 天 空中 点 源 位 置 与 6 RARER 
去 卷 积 的 残 差 图 组 成 . HOGBOM-CLEAN 算法 是 
一 个 很 好 的 近似 方法 , 该 算法 基于 点 源 模型 , 对 于 
展 源 的 去 卷 积 效果 不 佳 .MS-CLEAN 算法 对 此 提 
出 了 改进 , 该 算法 在 HOGBOM-CLEAN 算法 的 基 
MWE, 将 图 像 分 解 为 多 个 尺度 的 部 分 , 并 针对 每 个 
部 分 模拟 使 用 不 同 高 斯 源 的 PSF 进行 处 理 , 这 种 多 
尺度 的 方法 能 够 更 好 地 处 理 展 源 , 并 改善 了 PSF 29 
WORST] DU SA, 不 过 由 于 涉及 到 多 尺度 的 处 理 , MS- 
CLEAN 算法 需要 更 长 的 处 理 时 间 和 较 多 的 计算 资 
源 ， 虽然 CLEAN 算法 在 射电 天 文 观测 数据 处 理 
方面 取得 了 一 定 的 成 果 , 但 它们 无 法 从 根本 上 解决 
UV 频 域 不 完全 采样 的 问题 , 使 得 从 有 限 观测 数据 
重 构 可 信 的 图 像 成 为 一 大 挑战 . 新 一 代 射 电 望 远 
镜 , 如 低频 射电 望远镜 阵列 (Low Frequency Array, 
LOFAR) 9! 和 平方 公里 阵列 (Square Kilometre Ar- 
ray，SKA) O, 因 其 高 灵敏 度 、 高 分 辩 率 和 更 快 的 
这 天 速度 , 在 数据 处 理 方面 提出 了 新 的 挑战 ， 这 些 
射电 望远镜 产生 了 海量 的 观测 数据 四, 需要 快速 处 
J8. 为 了 应 对 这 些 挑战 , 迫切 需要 开发 新 的 算法 . 这 
些 算法 不 仅 要 能 处 理 海 量 的 数据 , 还 要 能 够 准确 地 
重建 图 像 , 并 能 对 成 像 过 程 中 的 不 确定 性 进行 量化 . 
近年 来 , 压缩 感知 技术 为 射电 干涉 成 像 领域 提供 了 
一 种 有 效 的 解决 方案 . 2008 年 Cornwell 等 人 在 论 
文中 提出 压缩 感知 技术 射电 成 像 的 可 行 性 外, 2009 
年 压缩 感知 应 用 于 射电 干涉 成 像 , 研究 结果 表明 该 
方法 的 可 行 性 及 其 相对 于 标准 干涉 成 像 技 术 的 优越 
性 BB,9， 在 此 基础 上 ,Dabbech 和 Oleg Smirnov 
等 人 结合 稀疏 平均 加 权 算 法 上 9， 提出 了 一 种 自 适 
应 预 处 理 原 始 对 偶 算 法 (Adaptive Preconditioned 
forward-backward Primal-Dual, PPD) !J, 该 算法 
的 目的 在 于 估计 VLA 观测 数据 中 校准 误差 和 未 知 
的 噪声 水 平实 验 结果 表明 , 与 MS-CLEAN "! 相 
比 , 自 适 应 PPD 算法 展现 了 更 优越 的 性 能 ，Li 等 
人 使 用 澳大利亚 平方 公里 阵列 探 路 者 (Australian 
SKA Pathfinder, ASKAP) 的 模拟 观测 数据 , 研究 
了 基于 健 里 叶 变 换 基 (Partial Fourier, PF) 和 各 向 
同性 非 抽 取 小 波 基 (Isotropic undecimated wavelet 
transform, IUWT) 的 去 卷 积 算法 1, 他 们 进行 了 
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两 种 不 同 稀 玻 基 下 展 源 图 像 重 构 的 质量 对 比分 析 ， 
并 在 IUWT 小 波 基 下 使 用 快速 迭代 阅 值 收缩 算法 
(A fast iterative shrinkage-thresholding algorithm, 
FISTA) "9! 将 图 像 复原 ,结果 表明 IUWT-FISTA 
算法 重建 的 射电 天 文 图 像 优 于 HOGBOM-CLEAN 
和 MS-CLEAN,FISTA AA ERER RE 6 1E 
则 化 问题 , 通过 梯度 下 降 的 结果 和 上 一 步 的 估计 进 
行 比 较 , 根据 步 长 和 残 差 的 差异 来 快速 迭代 更 新 估 
计 值 , 不 过 , FISTA 算法 对 输入 图 像 中 的 噪声 较为 
BUR, 各 输入 图 像 噪声 较 高 , 算法 可 能 无 法 有 效 去 
除 噪声 .为 了 处 理 新 一 代 的 射电 望远镜 阵列 , 更 精 
细 的 复原 射电 图 像 细 节 , 本 文采 用 了 原始 对 偶 分 裂 
算法 04 对 射电 干涉 成 像 的 逆 问 题 进行 改进 . 第 2 节 
介绍 了 处 理 观测 数据 的 稀 艳 优化 问题 , 包括 压缩 感 
知 框架 下 的 射电 干涉 成 像 逆 问 题 及 本 研究 所 采用 的 
图 像 复原 算法 ， 第 3 节 首 先 模 拟 了 SKAL-low 观测 
数据 ,然后 对 其 进行 了 实验 验证 ,对 比分 析 了 不 同 
算法 在 复原 效果 上 的 结果 , 最 后 对 本 算法 进行 深入 


讨论 和 总 结 . 


2 压缩 感知 框架 下 射电 阵列 干涉 成 像 逆 问 


H 


题 


2.1 稀 统 优化 问题 的 一 般 形 式 


Candes 等 人 的 压缩 感知 理论 研究 中 
提 h 15718], 在 一 定 条 件 下 ， cep 


xER” = {xh WDA xb — + Fp ot E Hu dE 
We RWN = (Pi) |, MARR ERAF V 
下 , 信号 x 可 以 分 解 为 ge RY = lada, 为 此 信 
号 x 可 以 定义 为 x= ya, 其 中 只 有 少量 有 意义 的 系 
数 oi, 大 部 分 系数 为 零 , 即 信号 x 含有 的 非 零 元 数 
H K « N, 考虑 背景 噪声 ne R" = {nhm 获得 
M 个 线性 观测 值 可 以 表示 为 : 


y=Oatn (1) 


yeR^z (Yih ckem 为 观测 值 , 其 中 O = Hy 是 一 个 
从 RY 到 RY 的 降 维 矩阵 , H e RON Rea TEAR Ri iC 
表示 下 的 降 维 矩阵 ， 从 公式 (恢复 a 需要 求解 以 


XXX-2 


EE XASEUE 
min llallos.tlly — Oall, < € (2) 
aeRN 


公式 (2) 是 一 个 非 确 定性 多 项 式 (Non-deterministic 
Polynomial, NP-hard) 问题 , 通常 使 用 €, 范 数 代 
Tk ly 范 数 将 其 转换 为 一 个 凸 优 化 问题 9, 将 其 等 
价 于 : 


minllall,s.tly — Oall, < € (3) 
aeRY 


2.22 ”射电 阵 干涉 成 像 逆 问题 


视 场 比较 小 时 ， 源 离 观测 相位 中 心 (B 7 = 
0,m = 0) 非常 近 , 根据 Van Cittert-Zernike 定理 ， 
可 见 度 函数 Vu, v) 和 天 空 亮度 分 布 1(1,m) 是 一 对 
二 维 健 里 叶 变 换 关系 Po, 


图 1 UV 平面 投影 基线 的 示意 图 


Figure 1 Schematic diagram of UV plane projection baseline 


如 图 1 所 示 ， 两 个 射电 望远镜 构成 一 条 物理 基 
线 ， 从 源 的 角度 看 向 物理 基线 可 以 得 到 投影 基线 ， 
物理 基线 的 位 置 和 天 空中 的 观测 方向 共同 决定 了 这 
条 投影 基线 的 特征 , 而 投影 基线 则 决定 了 所 测 天 空 
区 域 的 UV 空间 分 布 ， 随 着 地 球 自转 , 观测 时 间 
的 增加 , 在 UV 平面 形成 一 条 轨迹 . 每 两 个 射电 望 
远 镜 构 成 一 条 物理 基线 , 对 应 于 一 系列 可 见 度 函 数 
V(u, v) 的 采样 点 . 全 部 基线 所 获得 的 复 可 见 度 函数 
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与 亮度 分 布 之 间 有 以 下 关系 : 
V(u, v) = Í T I(l, m)e 7" didm (4) 
S 


对 公式 (4) 两 边 进行 传 里 叶 逆 变换 即 可 得 到 脏 图 
lul, m), 与 传统 的 奈 奎 斯 特 采 样 不 同 , 脏 图 是 “ 真 
实 信 号 ”与 观测 系统 PSF 卷 积 的 结果 , 存在 噪声 的 
情况 , 可 表示 为 : 


Tairty = PSF * (lirue + n) (5) 


Kp, * 表示 卷 积 符号 , 成 像 即 为 通过 采样 得 到 
的 已 知 可 见 度 数据 和 采样 函数 重建 未 知 的 原始 天 空 
图 像 1. 实际 观测 数据 受到 噪声 n. 误差 和 干扰 的 
影响 , 导致 图 像 模糊 , 为 了 获得 高 质量 的 成 像 结 果 ， 
需要 对 数据 进行 校准 和 去 卷 积 , 以 抑制 或 消除 脏 束 
SSNS. 在 射电 天 文 图 像 去 卷 积 问题 中 ,可 
以 求 出 很 多 个 射电 图 像 解 与 PSF 进行 卷 积 后 的 脏 
图 , 这 是 一 个 病态 问题 . 在 此 情况 下 , 可 以 通过 图 像 
的 先 验 知识 作为 限制 条 件 , 从 众多 的 解 中 选 出 最 合 
适 的 射电 图 像 解 . 因此 , 本 文 结合 最 新 的 压缩 感知 
理论 研究 一 种 新 的 解决 方法 . 假设 观测 目标 是 点 源 ， 
将 基于 压缩 感知 框架 的 重建 方法 应 用 于 此 , 从 而 公 
式 (4) 中 的 射电 天 文成 像 的 问题 可 以 描述 为 : 


min lh, stV = MinaseF 1 (6) 


如 图 2 所 示 , 矢量 工 表示 待 观测 目标 的 空间 亮度 图 ， 
矢量 V 为 通过 天 线 阵列 获得 的 观测 数据 , 由 于 UV 
平面 欠 采 样 , 其 中 包含 很 多 零 元 素 ; 矩阵 HAE 
叶 变 换 矩 阵 ; Mnasx 表示 二 进 制 掩 膜 矩 阵 , 其 中 1 表 
示 相 应 位 置 有 观测 数据 , 0 表示 没有 观测 数据 ， 针 
对 有 噪声 的 情况 , 该 问题 可 进一步 表述 为 : 


min ||/|I,, S.t ||Mmask FI = VII X 2 (7) 


其 中 , 6 表示 观测 数据 的 不 确定 性 , BID Sl PE 
对 于 展 源 目 标 , 上 述 方法 可 能 不 适用 . 通过 分 析 射 
HT REY Bil o TTI, 发 现 IUWT 小 波 字典 
是 一 个 有 效 的 射电 天 文 图 像 稀 朴 表示 的 图 像 模型 ， 
本 文 将 使 用 该 方法 进行 数学 建 模 . 其 中 , WON Ri it 
表示 变换 矩阵 , 其 逆 变换 为 W^, 基于 IUWT 小 波 
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< 
= 


mask 
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UV Coverage 


Measurement Matrix 
(Masking - Fourier) 


H =M mxf 


mask 


图 2 成 像 逆 问题 的 示意 区 
阵 , I 是 原始 天 空 模型 , 可 在 字典 WP TEN o. 


I =W} a 


, 其 中 Y 是 可 见 性 数据 , Minask 对 应 系统 采样 函数 , F. 表示 传 里 叶 变换 矩阵 , H 表示 是 降 维 变换 矩 


Figure 2 Schematic diagram of the imaging inverse problem,where V is the visibility data, Mmask corresponds to the system's 


sampling function, F denotes the Fourier transform matrix,H signifies a dimensionality reduction transformation matrix,/ is the 


original sky model, which can be sparsely represented in the dictionary W as « 


字典 和 压缩 感知 的 去 卷 积 方法 TUWT-CS 可 以 表示 
为 : 


minllalli, s.tV = Minask Wa (8) 


EH, o-WI 表示 的 是 经 WT 小 波 字典 表示 
AMARA, CAR UUM A re RE ESL RIA tu de 


23 图像 复原 算法 


在 射电 天 文 图 像 处 理 领域 , 小 波 变换 被 广泛 认 
为 是 一 种 非常 适用 的 方法 多， 这 种 方法 特别 适合 
处 理 表现 出 各 向 同性 或 准 各 向 同性 的 天 体 图 像 , 例 
如 恒星 、 星 系 或 星系 团 等 ，IUWT 小 波 变 换 进行 
到 像 处 理 的 具体 步 又 如 下 : 通过 对 一 幅 图 像 x 进 


m. M mask 表示 测量 矩阵 ， HS ri fern A E d] E 
矩阵 有 限 等 距 性 质 (Restricted Isometry Property, 
RIP). 针对 有 噪声 的 情况 , 有 : 


min llall, s.t| Ms FW"! — Vl < € (9) 
Tei o RUP IC TEARS TE P. YS CY bit 7G 
验 正则 化 约束 条 件 可 以 得 到 更 精确 的 解 22:291. 根据 
式 (9), MAGICA Aa, 问题 可 以 改写 成 拉 格 
W HJE: 


arg min awal], + | Mss FW ~'a -V 


i (10) 
从 贝 叶 斯 理论 的 角度 看 ， 式 (10) 中 4 可 视 
为 平衡 先 验 模型 Wah FR KR EROR 
IIM mass a — VIE 的 参数 . 


行 分 解 , 可 以 表示 为 x*=x/ + È w, 其 中 o! 表示 

j=1 
EREDE. wigo 表示 不 同 尺度 的 小 波 分 量 , 在 
这 里 , 第 一 级 小 波 分 量 (j=1) 对 应 于 最 高 频率 的 细 
节 部 分 ， 在 本 研究 中 , 为 了 凸显 图 像 局 部 细节 和 高 
频 特 征 , 选择 使 用 J—3 并 忽略 粗 尺度 分 量 , 基于 公 
式 (3) 和 公式 (10), 稀 琉 正则 化 的 最 小 化 问题 可 以 表 
述 为 : 


arg min 2 ly ~ Half + |w o x]. (11) 
其 中 x=ya, y 实现 了 没有 粗 尺 度 的 小 波 变换 ， 
Wk) 是 一 个 加 权 和 矩阵 ,& 值 范 围 为 0 和 ES kmax, © 
表示 Hadamard FR. W(k) All x PRE J x pti FE 
EE. k 是 重新 加 权 指 数 . €, 范 数 通常 通过 软 阅 值 实 
现 , 表示 为 kl, = bs, 软 阔 值 操作 符 具 体形 式 如 
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F: 


, Ix;l >A 


, 其 他 


STa(x;) = 


h — Asign(x;) (12) 


ERRER APPS BE II RN, 
而 包含 较为 重要 信息 的 较 大 小 波 系数 则 会 有 所 减 
小 ， 这 会 导致 信号 重建 这 程 中 出 现 一 定 的 偏差 , 为 
了 优化 这 种 不 平衡 并 更 接近 Co 范 数 ,可 通过 重新 
加 权 方法 来 改进 中, 记 $* 为 方程 的 解 , 对 式 (11) 中 
W*, 在 给 定 的 值 下 , 权重 的 递归 关系 定义 为 : 


1 


(k+1) _ (k) 
We = W (13) 


其 中 WO (XI pg. OAPI SR 
在 根据 去 卷 积 过 程 中 传播 到 估计 的 射电 图 像 不 确定 
性 来 调整 小 波 系 数 . 其 目的 是 为 受 观测 噪声 严重 影 
响 的 小 波 系数 分 配 较 大 的 权重 , 为 其 他 小 波 系数 分 
配 较 小 的 权重 , 同时 , 通过 重新 加 权 , 减少 了 对 那 
些 远 高 于 其 预期 噪声 水 平 小 波 系 数 的 惩罚 , 从 而 减 
轻 了 由 6 范 数 引起 的 偏差 .对 于 小 波 和 矩阵 ， 可 以 
表示 为 y= (WT, us EP yia 表示 图 像 x 在 
第 j 个 小 波 尺度 下 未 经 抽样 的 小 波 变 换 ， 当 PSF 
卷 积 和 矩阵 H nput A AN PY) RIP 条 件 
时 , 可 以 通过 对 向 量 Hy 进行 羡 值 处 理 来 直接 实 
现 去 噪 和 复原 图 像 ， 为 此 , 在 尺度 j FEX NE 
向 量 为 # = kolo" |, 其 中 o 是 射电 图 像 中 噪 
声 的 标准 偏差 , k 是 与 尺度 相关 的 调节 参数 ， 在 
噪声 服从 高 斯 分 布 的 情况 下 , k; 可 以 取 3 或 4 以 
获得 更 佳 结果 .为 了 在 稀 下 域内 产生 更 加 稀 芯 的 
解 ， KLARA TREERE AR) , 同时, BÀ 
值 A 的 选择 将 根据 问题 类 型 的 不 同 而 变化 , 在 去 
噪 情 况 下 ,可 以 设置 4 = noew, 其 中 es Im 
中 噪声 标准 偏差 的 估计 值 . 假设 整个 场 遵 循 高 斯 
白 噪 声 分 布 , 可 以 通过 以 下 公式 估计 噪声 标准 偏差 : 
Tes = (1.4826) x MAD(Y), 其 中 MAD 表示 中 值 绝 
对 偏差 ， 定 义 为 MAD = median(|x; — median(x;))). 
n 是 一 个 乘法 因子 , 对 应 于 # 中 的 k; Alo. 根据 
公式 (11) 和 公式 (13)， 可 以 将 问题 表述 为 以 下 优化 


1) https://www.python.org/ 
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问题 


ea 


(14) 


直接 求解 公式 (14) 是 不 现实 的 , 正如 图 2 所 表示 的 卷 
积 是 奇异 矩阵 ,没有 精确 的 逆 矩 阵 ， 这 将 无 法 对 去 
卷 积 计算 ， 为 此 ，Condat 算法 通过 在 目标 函数 中 
引入 正则 化 项 (惩罚 项 ), 使 得 在 求解 拟 合 模型 时 引 
入 额外 信息 ,以 防止 过 拟 合 和 观测 和 矩阵 非 满 秩 的 情 
况 , 可 以 用 Hz 来 近似 H, 得 到 # = kei |. 
从 而 初始 权重 将 定义 为 WO = [TT OWT 
解决 上 述 讨论 中 提出 的 最 优化 问题 ， 本 文采 用 了 
Python" gj fi& 53€ 3:34] Condat 在 其 文中 3.1 部 
分 描述 的 原始 对 侦 分 裂 算法 5. 在 该 算法 中 , 忽略 
了 误差 项 , 其 原始 凸 优化 问题 可 以 描述 如 下 : 


arg min F(x) + G(x) + K(L(x)) (15) 


其 中 函数 F 是 一 个 梯度 为 VE 的 凸 函数 , 而 G 和 
K 则 是 可 以 通过 封闭 形式 有 效 求解 的 近似 算 子 . L 
是 一 个 计算 其 正则 化 和 伴随 操作 的 有 界线 性 算 子 . 


算法 1 选择 近 端 参数 r > 0,0 > 0、 正 松弛 参数 pn 和 初始 估计 值 (xo, Yo). 
然后 对 每 个 上 > 0 进行 迭代 


Žn+1 = proxig(x, — TVF (Xn) — TL'(y4)) 
na = proxgek- On 十 TL(2Xn+41 = Xn)) 
(Xn; Yn+1) = pina aa) t (1 — Pn) (Xn. Yn) 


在 算法 1 中 , FOAL 了 分 别 表示 原始 变量 和 对 
(BE, 其 中 prox O) =y- oproxe (2). 4S 
收敛 后 , 原始 变量 将 对 应 于 最 终 解 ， 即 去 卷 积 后 的 
射电 图 像 ， 如 图 3 所 示 , 该 算法 的 实现 需要 输入 图 
像 数 据 、PSF、 适 代 次 数 n 和 重新 加 权 次 数 me 
其 中 原始 近 端 算 子 proxrc 表示 正则 化 项 的 正 性 约 
束 ， 确 保重 建 图 像 的 像素 值 保持 非 负 ， 梯度 算 子 
VF (a) = &'H" (Hoa - y) 在 每 次 迁 代 中 计算 当前 重 
建 图 像 相 对 于 损失 函数 的 梯度 ， 这 个 梯度 会 被 用 
于 调整 重建 图 像 ， 使 其 在 梯度 方向 上 减 小 损失 函 


Fi 
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数 的 值 . 对 偶 近 端 算 子 proxa(ž) 对 应 于 方程 中 
权重 Wk) 的 软 阔 值 稀疏 处 理 ， 线 性 运算 符 L 对 
应 于 小 波 变换 B， 在 式 (14) 中 , 8 中 的 k; 设置 为 
ki = 3,k2 =3,k; = 4, 默认 重新 加 权 次 数 设 置 为 1. 
对 侦 近 端 算 子 参数 设置 为 T=0=0.5. 对 于 所 有 实现 ， 
松弛 参数 设置 为 p,=0.8. 计算 成 本 函数 , BERS 
间 成 本 函数 的 变化 小 于 0.0001 时 , 算法 被 认为 收 
SX. 当 迭 代 次 数 和 重 加 权 次 数 满足 条 件 时 , 输出 结 
果 . 


输入 脏 图 和 PSF, 设置 迭代 
RA n HERRERAN, . 


RERA X FOU y» UM, 定义 梯度 算 子 Y 和 线 


RSL, REMAS prox c 和 对 偶 近 端 算 子 


prox, 分 别 为 正则 项 和 软 阔 值 项 . 


结合 梯度 算 子 、 近 端 算 子 、 对 偶 近 端 算 子 V. 


prox. prox, 计算 成 本 函数 ( tolerance x 107* y 


结合 x. y. V. L, prox. prox, AMAR 


数 计算 condat 函数 (7 = 0.5,0 = 0.5, p, = 0.8). 


迭代 和 更 新 权重 值 
达到 n 和 ne? 


图 3 重 构 天 空 模型 的 算法 流程 
Figure 3 Workflow of algorithm for reconstructing the sky 


model 


2) https://casa.nrao.edu/ 

3) https://github.com/OxfordSKA/OSKAR 

4) https://casaguides.nrao.edu 

5) https://code.google.com/archive/p/csra/downloads 
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3 实验 验证 


3.1 ”测试 方法 与 结果 


本 文 使 用 通用 的 天 文学 软件 包 ”(Common As- 
tronomy Software Applications, CASA) PS, 4% 
大 学 开发 的 SKA 模拟 软件 包 (Oxford SKA Radio 
Telescope Simulator, OSKAR)? fil Python 编程 语 
言 来 实现 上 述 算 法 .为 了 评估 算法 的 去 卷 积 性 能 
并 排除 其 他 因素 对 测试 结果 的 影响 , 制作 了 原始 天 
空 模型 , 主要 模拟 了 两 个 点 源 和 位 于 猫 大 座 的 M51 
涡 状 星 系 . M51 原始 的 FITS 参考 图 像 数 据 可 在 
NRAO HI FARO. 天 空 模型 的 参数 设置 与 参考 图 
像 生成 相对 应 , 包括 赤道 坐标 中 的 位 置 ( 赤 经 和 赤 
纬 )、Stokes 参数 (I. Q, U., V), 在 参考 频率 下 工 的 
流量 密度 、 谱 指数 、 旋 转 测量 以 及 其 他 天 空 模型 参 
数 默 认 设置 为 零 ， 对 于 测量 集 数 据 (Measurement 
Set, MS) 格式 的 天 文 数据 , HOGBOM-CLEAN 和 
MS-CLEAN 算法 适用 , 并 使 用 CASA 软件 包 进 行 
分 析 : 而 IUWT-FISTA 和 IUWT-CS 算法 则 专门 
针对 FITS 格式 的 天 文 数据 其 中 IUWT-FISTA 
使 用 Matlab m5, IUWT-CS 使 用 Python 实 
W. 在 本 文中 , 模拟 测 站 为 射电 望远镜 阵列 模型 
SKA1l-low L271， 中 心 参考 位 置 为 [116.764, -26.825], 
如 图 4 所 示 , SKAl-low 共计 131072 个 对 数 周期 天 
线 , 分 成 512 个 阵列 站 组 成 , 每 个 站 包含 256 个 天 
线 , 其 中 296 个 位 于 到 中 心 核心 区 域 , 其 他 216 个 沿 
着 3 个 旋 臂 部 署 . 通过 设置 M51 在 中 心 相 位 位 置 
RA=20°, Dec=—30° 处 , 使 用 6 个 通道 (100、110、 
120、130、140 和 150MHz) 进行 观测 模拟 , 参考 频 
AM. 100MHz 开始 , 以 10MHz 的 带宽 递增 ,观测 
数据 的 数值 属性 详 见 表 1. 

观测 过 程 包括 4、6、8、12 小 时 的 持续 时 间 , 分 
别 进行 了 24 个 快照 . 不 同 观 测 时 长 的 干涉 模拟 耗 
时 分 别 为 22413.36、22623.58、22352.54、20819.42 
秒 ， 通 过 应 用 程序 oskar_imager， 对 生成 的 可 见 
度 文 件 进行 处 理 ， 得 到 的 图 像 覆 盖 了 观测 相位 中 
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表 1 模拟 干涉 数据 参数 设置 


Table 1 Parameter Settings of Simulate Interferometric Data 
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Parameter Settings M51 Two point sources M31 

Start Frequency (Hz) 1.0e+08 1.5e+08 1.0e+08 

Number Of Frequency channels 6 1 3 
Frequency Increment (Hz) 1.0e+07 1.0e+07 1.0e+07 
Start time (UTC) 15-08-2023 12:00:00 15-08-2023 12:00:00 15-08-2023 12:00:00 

Observation length 4h/6h/8h/12h 12h 12h 

Number of time steps 24 24 24 

Channel Bandwidth (Hz) 300 300 300 

Time Average 10 10 10 

Output Data Format MS/FITS MS/FITS MS/FITS 


心 4° MA, UV 覆盖 情况 如 图 5 所 示 ， 左 侧 两 列 
观测 开始 时 间 为 12:00:00, 最 长 基线 可 达 65850m, 
右 侧 两 列 为 局 部 UV 覆盖 图 (为 了 选用 覆盖 率 
更 高 的 模拟 观测 数据 ， 本 文选 择 12h 的 数据 进 
行 算法 测试 ， 所 得 到 的 脏 图 如 图 6 所 示 ， 图 像 大 小 
为 256*256 像素 ,单位 为 Jy/beam, 频率 通道 为 
1~6( 即 100MHz-150MHz). 图 像 亮度 值 范围 在 [0， 
1.25] 区 间 ， 最 大 流量 为 1.3184Jy， 通 过 应 用 程序 
oskar vis to ms, 将 可 见 度 文件 转换 为 MS 格式 ， 
以 便 在 CASA 软件 包 中 进行 后 续 处 理 . 


SKA1-low Full Stations 


SKA1-low Core Area 


North (km) 
North (m) 


~~ East (km) 


SKA1-low Core Stations 


North (m) 


~~ East (m) 


East (m) 


图 4 SKA1l-low 阵列 布局 图 
Figure 4 Schematic of SKA1-low array layout. 


当 考 虑 点 源 重 建 时 ，CLEAN 算法 通常 是 首选 


方法 . 


因此 , 本 研究 首先 比较 了 IUWT-CS 方法 与 


HOGBOM-CLEAN 算法 在 处 理 这 类 简单 源 时 的 性 
能 . 通过 模拟 了 两 个 1Jy 的 点 源 , 采用 SKA1-low 
进行 了 12 小 时 的 观测 (参见 表 1 中 的 数据 集 ) . 两 个 
点 源 中 , 一 个 位 于 相位 中 心 , 男 一 个 位 于 距离 相位 
中 心 5” 到 30” 之 间 . 这 样 设置 可 以 比较 PSF 卷 积 
点 源 后 的 分 离 程度 , HOGBOM-CLEAN 和 IUWT- 
CS 性 能 的 比较 分 为 三 个 区 域 : 源 未 解析 、 源 部 分 解 
析 、 明 确 解 析 . 由 于 源 靠近 相位 中 心 , 第 2.2 节 开始 
提 到 的 效应 可 以 忽略 .在 观测 数据 集中 , 通过 加 入 
了 不 同 频率 (1.49e8, 1.50e8, 1.51e8) 对 应 不 同 的 
噪声 水 平 (0.1~0.3 高 斯 随机 和 白 噪 声 ) , 大约 达 到 了 
10 的 信 噪 比 (S/N) 水 平 , 通过 使 用 CASA 包 中 的 
HOGBOM-CLEAN 和 Python 实现 的 IUWT-CS 
方法 进行 了 去 卷 积 . 对 于 CLEAN 算法 , 通过 执行 
T 20000 次 迭代 , 并 设置 briggs=0.5, imsize=256. 
IUWT-CS 执行 了 200 次 迭代 ， 更 新 权重 次 数 1 
次 , 图 7 展示 了 获得 的 结果 : HOGBOM-CLEAN 重 
建 ( 左 栏 )、IUWT-CS 重建 (PE) 以 及 脏 图 (A 
RD). 各行 对 应 9 值 的 不 同 , HOGBOM-CLEAN 
的 天 空 模型 通常 由 许多 点 pixel 像素 组 成 ， 有 两 
个 最 大 值 的 点 源 像素 值 , 两 个 源 完全 可 分 辩 开 来 变 
成 大 约 2~3 个 点 pixel 像素 组 成 , 与 HOGBOM- 
CLEAN 不 同 , IUWT-CS 理论 的 天 空 模型 即 为 两 
MEA 1 的 点 pixel， 因 为 在 有 限 的 分 辨 率 下 , E 
缩 感 知 框架 下 的 输出 结果 即 为 真实 天 空 的 最 佳 表 
示 [ 妆 ， 模 拟 观测 的 脏 束 尺寸 约 为 15.6”x9.7”， 当 
A6—5" Ht, HOGBOM-CLEAN 和 IUWT-CS 并 


XXX-7 


张 讯 等 . 


(a) x10 UV Coverage(4h) (b) x10 UV Coverage(6h) (e) x10 UV Coverage(4h) UV Coverage(6h) 
6 6 i 
4 4 
2 
2 2 
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> > > 
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-4 -4 
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-6 -5 
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(c) «104 UV Coverage(8h) (d ) x10* UV Coverage(12h) UV Coverage(8h) x10? UV Coverage(12h) 
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6 
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E -2 
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-6 -h 
-6 -4 -2 0 2 4 6 -6 -4 -2 0 2 6 4 2 4 
| U Imi x10% U (m) x10* Ulm) x10? Ulm) x10? 
图 5 连续 观测 4. 6. 8. 12 小 时 时 二 维 频率 空间 UV 平面 覆盖 , 图 (a). K (b). 图 (c) 和 图 (d) 开始 观测 时 间 为 12:00:00， 


x 轴 取 值 区 间 为 (-71254.425m, 71254.425m), y 轴 取 值 区间 为 (-68238.775m, 68238.775m). 图 (e), 图 (£). 图 (g) 和 图 (h) 


开始 观测 时 间 为 12:00:00, x 轴 取 什 


区 间 为 (-500.00m, 500.00m), y 轴 取 值 区 间 为 (-500m, 500.00m). 


Figure 5 The UV-plane coverage in two-dimensional frequency space for continuous observations of 4, 6, 8, and 12 hours, respec- 
tively. Figures (a), (b), (c), and (d) start observation at 12:00:00, with the x-axis ranging from -71254.425m to 71254.425m and the 
y-axis from -68238.775m to 68238.775m. Figures (e), (f), (g), and (h) commence observation at 12:00:00, with the x-axis spanning 
from -500.00m to 500.00m, and the y-axis ranging from -500m to 500.00m. 


J2000 Declination (degrees) 
—31°00'00" —30*00'00" —29*00'00" —28*00'00" —31*00'00" —30*00'00" —29*00'00" —28°00'00” 


130MHz, 140MHz, 150MHz 


01^24"005 01^20"00* 01^167"00* 01^12700* 01^24"005 01^20"00* 01^16"005 01^127"00* 


01^24"005 01^20"00* 01^16"00* 01^12700* 


J2000 Right Ascension (hours) 


图 6 观测 时 长 为 12h 时 , 6 个 不 同 频率 通道 上 观测 数据 生成 的 脏 图 , 从 左 到 右 从 上 到 下 分 别 为 100MHz , 110MHz, 120MHz, 


Figure 6 Dirty maps generate from observational data across six different frequency channels during 12-hour observation period, 
arranged from left to right and top to bottom for frequencies 100MHz, 110MHz, 120MHz, 130MHz, 140MHz, and 150MHz respec- 


tively. 
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图 7 SKAl-low 观测 点 源 及 重建 结果 , HOGBOM-CLEAN( 左 列 )、IUWT-CS( 中 间 ) 和 脏 图 ( 右 列 ). 
Figure 7 Observational and reconstruction results of point sources using SKA1-Low, featuring HOGBOM-CLEAN (left column), 


IUWT-CS (middle), and dirty maps (right column). 

未 区 别 出 两 个 源 ， 对 源 的 全 宽 半 高 (Full width at 
half maximum, FWHM) 进行 椭圆 拟 合 给 出 了 源 的 
大 小 , 在 这 种 情况 下 , HOGBOM-CLEAN 去 卷 积 
后 源 的 尺寸 为 16”x10”, IUWT-CS 为 12”x8”, 4 
A6—15" Hf, 这 两 个 源 的 位 置 不 确定 性 仍然 较 高 , 但 


高 S/N 条 件 下 , 范围 在 [10, 100], IUWT-CS 的 分 
辨 性 超过 HOGBOM-CLEAN 的 2-3 fi. 24 S/N 
降低 时 , IUWT-CS 的 可 分 辨 性 趋向 于 HOGBOM- 
CLEAN 的 表现 . 除了 角 分 辩 率 误差 外 , 通过 查看 重 
建 源 的 通 量 密度 误差 受到 背景 噪声 水 平 的 影响 . 对 


当 A0 增加 到 30” 时 , 两 个 源 之 间 的 测量 和 流量 密 
度 误差 开始 减 小 , HOGBOM-CLEAN 7 15"x10", 
IUWT-CS 方法 为 T"x4". 且 在 源 分 离 方面 的 效果 
优 于 HOGBOM-CLEAN 算法 ， 随 着 Ao 的 增加 ， 
最 后 两 个 源 在 正确 的 位 置 被 明确 分 开 ， 在 给 定 适 
当 的 Ae 和 信 噪 比 水 平 下 , 通过 从 不 同 的 信 品 比 水 
平 下 分 辨 两 个 不 同 源 的 有 效 分 辨 率 ， 图 8 为 通过 俩 
用 HOGBOM-CLEAN 和 IUWT-CS 算法 获得 的 
结果 曲线 .对 于 每 个 S/N, 记录 两 个 源 分 离 的 角度 
A6 来 证 明 IUWT-CS 带 来 的 有 效 改进 . 结果 表明 ， 
HOGBOM-CLEAN 的 可 分 辨 性 对 S/N 的 依赖 性 
AR, 这 使 其 成 为 在 各 种 S/N 条 件 下 检测 点 源 仅 是 
一 个 稳定 的 复原 算法 .IUWT-CS 的 表现 不 同 , 在 


T IUWT-CS, 它 在 低 噪声 条 件 下 为 2%, 在 高 噪声 
条 件 下 可 达 20%, 而 对 于 HOGBOM-CLEAN, 这 个 
范围 是 2% 到 1696. 总 的 来 说 , 点 源 情况 下 , IUWT- 
CS 复原 效果 几乎 与 HOGBOM-CLEAN 的 结果 一 
样 好 , 并 且 在 较 高 S/N 数据 下 提供 了 改进 的 角 分 辨 
率 , 这 表明 了 从 采样 不 足 的 干涉 数据 中 极 大 提高 图 
像 分 辩 率 的 可 能 性 . 
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=@ IUWT-CS 
=@- HOGBOM-CLEAN 


Angular separation (") 


图 8 重建 点 源 的 最 小 角 分 辨 率 值 


Figure 8 Minimum angular resolution values for the recon- 


struction of point sources. 


在 考虑 展 源 重 建 时 , 需要 注意 展 源 与 点 源 的 不 
I] AREE. 展 源 在 健 里 叶 变 换 和 矩阵 前 并 不 是 稀 蚊 
表示 的 , 因此 需要 对 其 进行 小 波 变换 , 此 外 , 使 用 不 
同 的 加 权 方 式 会 影响 点 扩散 函数 的 大 小 , 为 了 满足 
空间 分 辨 率 和 表面 亮度 灵敏 度 的 基本 要 求 , 并 生成 
适当 大 小 的 PSF, 本 次 测试 中 使 用 Briggs 加 权 方 
RPO, ERE 12 小 时 的 观测 时 间 在 100MHz 下 的 
脏 图 ，HOGBOM-CLEAN 的 参数 设置 如 下 : 网 格 
化 gridder='standard、 权 重 briggs=0.5、 迭 代 次 数 
niter 一 25000、 增 益 gain=0.1. 通过 实测 脏 图 的 噪声 
背景 区 域 , 发 现 最 低 均 方 根 值 约 为 0.06Jy/beam, 据 
此 计算 得 到 的 等 效 阐 值 threshold=60mJy /beam. 
从 图 9 第 一 列 可 以 看 出 , 去 卷 积 之 后 , 该 图 像 仍 然 存 
TEV HS 图 像 较为 模糊 .与 脏 图 6 相 比 , 图 像 结 
构 逐 渐 清 晰 ， 由 于 点 源 假设 模型 限制 ,其 天 空 模型 
MDA 6 函数 为 源 ， 无 法 本 质 上 逼近 扩展 结构 ， 从 
残 差 图 中 可 以 观察 到 源 的 流量 和 结构 并 没有 完全 被 
清除 , 误差 图 像 中 大 部 分 的 源 结构 依然 保留 ，MS- 
CLEAN 的 参数 设置 为 :scales=|[0 3 10 30]、 小 尺度 
Ani smallscalebias=0.9、 和 迭代 次 数 niter=25000, 
增益 gain=0.1、 阔 值 threshold=60 mJy/beam、 权 
重 briggs=0.5. 从 图 9 第 二 列 可 以 看 出 , MS-CLEAN 
相 比 HOGBOM-CLEAN, 图 像 更 加 清晰 , 残 差 医 
中 的 噪声 和 源 得 到 了 更 好 的 分 离 。 然而, 仍然 存在 
微弱 的 扩展 结构 MS-CLEAN 通过 将 原始 天 空 模 
型 表示 为 不 同 尺度 的 6 函数 组 合 , 改善 了 由 于 单一 


张 讯 等 . 


6 函数 模型 的 限制 但 误差 图 像 表明 , MS-CLEAN 
重建 了 大 部 分 的 扩展 结构 , 保留 了 一 部 分 微弱 结构 
的 扩展 源 , MS-CLEAN 只 将 图 像 分 解 成 更 多 尺度 
函数 来 复原 , 无 法 本 质 上 处 理 更 多 的 尺度 函数 来 表 
示 , 在 重建 扩展 结构 方面 还 存在 局 限 ， 而 IUWT- 
FISTA 和 IUWT-CS 的 尺度 函数 模型 能 够 有 效 解 
决 这 一 问题 . IUWT-FISTA 的 参数 设置 为 : 正则 化 
参数 lambda-0.00001, 3*fC6 XX niter=50, posi- 
tiveflag-1, waveletflg-1. /JWJEZK3E level=4. 观察 
到 9 第 三 列 , 相 比 于 图 6 中 的 脏 图 ， 去 卷 积 后 的 图 像 
有 明显 改善 , 没有 受到 旁 欠 的 影响 ,背景 噪声 略 有 
增加 . 残 图 中 大 部 分 源 基本 淹没 在 噪声 中 并 从 脏 图 
中 有 效 提取 出 来 , 相 较 于 天 空 模 型 ,去 卷 积 的 结 
中 局 部 的 扩展 源 变 得 更 加 清晰 平滑 . 这 由 于 FISTA 
算法 在 处 理 细节 部 分 时 表现 不 佳 , 将 噪声 误 判 为 图 
fea, 导致 结果 质量 下 降 ， 其 在 每 次 迭代 中 , 首 
先 计算 当前 估计 的 梯度 信息 ,然后 进行 阔 值 化 处 理 
以 获得 估计 的 图 像 , 最 后 通过 线性 插值 产生 新 的 佑 
计 结 果 . 虽然 这 能 够 快速 收敛 恢复 图 像 , 但 在 更 精 
细 控 制 估计 结果 和 受 噪声 影响 方面 不 如 Condat 处 
JH. IUWT-CS 的 参数 设置 为 : ERM niter—100, 
更 新 权重 迭代 weightniter==10、 小 波 水 平 level=4. 
迭代 结束 时 即 得 到 收敛 结果 . 更 新 权重 的 目的 是 在 
A (EK ARE Pro A EE, 确保 最 小 的 小 波 系 数 
不 会 被 置 零 ， 而 是 分 配 较 大 的 权重 .观察 图 9 第 四 
列 , 重建 结果 更 加 清晰 , 整体 上 去 除了 旁 激 的 影响 ， 
效果 明显 优 于 脏 图 , 结果 仅 有 微弱 的 模糊 ， 与 天 空 
模型 相 比 , 去 卷 积 的 结果 接近 真实 的 天 空 模型 图 像 . 
即使 在 源 较 弱 的 图 像 边缘 , IUWT-CS 的 结果 与 真 
实 图 像 非常 吻合 . 残 差 图 中 的 噪声 和 源 已 经 基本 分 
离 , 噪声 水 平 很 低 , 流量 和 峰值 剩余 流量 值 与 噪声 
水 平 相 差 不 大 ， 误 差 图 像 中 仅 剩 微弱 的 扩展 结构 . 
最 后 一 行 最 后 一 列 的 图 像 结果 (图 9), 数值 几乎 接 
近 零 值 . 


上 述 算法 均 在 配置 为 Ubuntu22.04 LTS 系统 
(12th Gen Intel(R) Core(TM) i9-12900H 2.50 GHz, 
处 理 核 数 16, 16GB 内 存 ) 上 运行 . 如 表 3 所 示 , 处 理 
一 张 256*256 像素 的 图 像 时 , 各 算法 的 运行 时 间 如 
F: HOGBOM-CLEAN, IUWT-FISTA 和 IUWT- 
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图 9 使 用 HOGBOM-CLEAN( 左 列 )、MS-CLEAN( 中 间 左 列 )、IUWT-FISTA( 中 间 右 列 ) fl IUWT-CS( 右 列 )， 从 模拟 
SKAl-low 阵列 观测 及 重建 M51 图 像 ， 从 上 到 下 : 洁 图 (Jy/beam)、 天 空 模型 (Jy/pixel)、 残 差 和 误差 图 像 . IUWT-CS 
恢复 图 像 和 模型 图 像 仅 在 缩放 方面 有 所 不 同 (前 者 使 用 HOGBOM-CLEAN 光束 的 光束 区 域 以 Jy/beam 为 单位 , 后 者 以 
Jy/pixel 为 单位 ). 洁 图 和 天 空 模型 数值 线性 归 一 化 为 最 大 流量 1, 误差 图 像 数 值 线 性 归 一 化 为 0.5, 残 图 图 像 的 色 阶 对 于 每 
一 列 都 是 不 同 的 , 如 右 侧 的 颜色 条 所 示 . IUWT-CS 重建 包含 从 数据 集中 恢复 的 高 空间 频率 信息 , 重建 结果 比 另外 三 种 算法 
获得 的 更 少 的 错误 结构 和 更 低 的 残余 水 平 . 

Figure 9 Simulate M51 image reconstructed from simulated SKA1-low observations using HOGBOM-CLEAN (left column), MS- 
CLEAN (middle left column), IUWT-FISTA (middle right column) and IUWT-CS (right column).From top to bottom: clean 
image (Jy/beam),sky model (Jy/pixel),residuals and error images. The IUWT-CS restored image and the model image differ only 
in scaling (the former uses the beam area of the HOGBOM-CLEAN beam in Jy/beam, the latter in Jy/pixel). The values of the 
clean image and the sky model are linearly normalized to a maximum flow of 1, and the values of the error image are linearly 


normalized to 0.5. The color scale of the residual image is different for each column, as shown in the color bar on the right. The 
IUWT-CS reconstruction contains high spatial frequency information recovered from the dataset, and the reconstruction results 
have less error structures and lower residual levels than those obtained by the other three algorithms. 

表 2 四 种 方法 重建 结果 的 数值 比较 


Table 2 Numerical comparison of reconstruction results using the four methods 


Algorithm HOGBOM-CLEAN MS-CLEAN IUWT-FISTA IUWT-CS 


RMSE 3.3770e-01 5.8880e-02 1.8325e-03 1.7466e-03 
SNR 5.1893 8.2820 1.4821 12.2219 
DR 2.9611 16.9836 545.6904 572.5319 
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表 3 四 种 方法 重建 不 同 尺寸 图 像 的 速度 比较 
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Table 3 Speed comparison for reconstructing images of different sizes using the four methods 


Image size(pixels) HOGBOM-CLEAN(s) MS-CLEAN(s) IUWT-FISTA(s) IUWT-CS(s) 
64*64 1.709 45.329 1.126 2.186 
128*128 1.717 77.463 2.290 3.569 
256*256 2.560 100.359 10.037 26.261 
512*512 5.963 124.908 63.112 68.108 
768*768 11.256 265.701 142.579 162.257 
896*896 17.292 407.929 220.130 264.825 
1024*1024 21.891 558.653 376.005 374.086 
ls 表示 算法 耗 时 秒 数 
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图 10 HOGBOM-CLEAN(Z: E). MS-CLEAN( E). IUWT-FISTA(ZE F) 和 IUWT-CSCR F), 四 种 不 同 重建 图 像 (y 


4H) 与 真实 图 像 (x 轴 ) 的 线性 关系 图 


Figure 10 Linear relationship graphs between the reconstructed images (y-axis) and the real images 


(x-axis) using four different 


methods: HOGBOM-CLEAN (top left), MS-CLEAN (top right), IUWT-FISTA (bottom left), and IUWT-CS (bottom right) 


CS 分 别 需要 2.560 fb. 10.037 PP. 26.261 fb. 而 
MS-CLEAN 则 需要 更 长 的 时 间 , 达到 100.359 fb. 
在 保持 上 述 环境 和 算法 参数 设置 不 变 的 情况 下 , 随 
着 图 像 大 小 的 改变 , 各 算法 的 处 理 时 间 也 会 有 所 变 
化 . 除了 HOGBOM-CLEAN, 其 他 算法 在 像素 大 
小 达到 512*512 之 后 开始 出 现 瓶 颈 , 消耗 时 间 显 著 
增加 ， 为 了 评估 重建 质量 , 使 用 信 噪 比 作为 重建 质 


量 指标 ， 


SNR = 201og i( — 


) (16) 


O xt 


其 中 oy 和 ous 表示 原始 天 空 图 像 的 标准 偏差 和 
误差 图 像 的 标准 偏差 之 比 go， 使 用 动态 范围 度量 


(Dynamic Range, DR), 


max (restored image) 
RMSE 


DR 是 指 重 建 图 像 的 峰值 流量 跟 残 差 图 像 标 准 差 
RMSE 之 比 H. 其 中 残 差 图 像 标准 差 定义 为 : 


2 
RMSE = >» (residual image) 
number of pixels 


DR = (17) 


(18) 


从 表 2 可 以 得 出 ， 相 对 于 HOGBOM-CLEAN 
和 MS-CLEAN 算法 , IUWT-FISTA 和 IUWT-CS 
算法 的 残 差 图 像 标准 差 小 了 一 个 数量 级 , 且 它们 的 
动态 范围 明显 更 高 ， 这 一 结果 表明 , IUWT-FISTA 
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和 IUWT-CS 算法 在 分 辨 较 弱 信和 号 源 方 面 有 优势 ， 
具体 来 说 , IUWT-CS 的 信 噪 比比 IUWT-FISTA 高 
出 了 4.74dB. 相对 于 MS-CLEAN, IUWT-FISTA 
的 信 噪 比 水 平 大 致 相同 ， 差 值 仅 为 0.79dB. 然而 ， 
HOGBOM-CLEAN 算法 的 信 噪 比 明 显 不 如 其 他 
算法 ， 在 所 有 算法 中 ，IUWT-CS 的 信 噪 比 最 高 ， 
比 HOGBOM-CLEAN 高 出 了 7.03dB, 比 IUWT- 
FISTA 高 出 了 4.74dB， 比 MS-CLEAN 高 出 约 
3.94dB. 这 些 数据 表明 ，IUWT-CS 在 去 噪 和 复原 
效果 上 相 比 其 他 算法 有 明显 优势 .图 9 展示 了 这 些 
算法 的 测试 结果 ,可 以 看 到 基于 IUWT-CS 的 重 
建 结 果 优 于 HOGBOM-CLEAN、MS-CLEAN 和 
IUWT-FISTA. 与 其 他 算法 相 比 , 基于 IUWT-CS 
的 洁 图 可 以 更 清晰 地 呈现 出 接近 真实 天 空 图 像 的 效 
R. 其 残 差 图 像 也 是 图 9 第 3 行 中 最 均匀 的 残 差 图 
ff, 表明 其 重建 图 像 相 比 其 他 算法 具有 更 少 的 伪 影 
以 及 更 广 的 动态 范围 . 在 理想 情况 下 , 真实 天 空 
素 与 重建 像素 之 间 的 关系 应 该 是 一 条 完美 的 线性 斜 
率 ， 只 对 受 噪 声 影响 的 像素 变 得 更 加 分 散 ， 图 10 中 
显示 了 全 分 辨 率 256*256 e AER E 其 中 
比较 了 重建 图 像 与 真实 天 空 图 像 之 间 的 关系 . 图 中 
的 y 轴 表 示 重 建 图 像 数 据 , x 轴 表 示 真 实 天 空 图 像 . 
从 图 中 可 以 看 出 , 基于 IUWT-FISTA fil IUWT-CS 
复原 算法 的 结果 更 接近 真实 的 天 空 图 像 . 
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根据 上 述 结果 分 析 , IUWT-CS 在 处 理 展 源 数 
据 时 表现 出 了 良好 的 复原 和 去 品 效 果 ， 然而 ,在 
实际 应 用 中 , 为 了 满足 不 同 的 科学 目标 和 任务 需求 ， 
需要 对 参数 进行 调节 以 获得 预期 的 图 像 质量 . 对 于 
展 源 , 选 定 IUWT 小 波 ， AVR BUA. A08 
数 、 近 端 原始 对 偶 算 子 等 参数 为 合适 的 取 值 情况 下 ， 
主要 考虑 有 两 个 关键 参数 会 影响 成 图 质量 : 一 是 迭 
代 次 数 , 二 是 更 新 权重 次 数 . 通过 调整 这 些 参 数 , 最 
终 复原 图 像 的 均 方 根 误差 、 动 态 范围 和 信 品 比 均 会 
发 生变 化 . 为 了 获得 图 像 的 最 大 动态 范围 或 信 噪 比 ， 
模拟 对 仙女 座 大 星系 M31 进行 观测 , 参数 设置 见 
Kl. 
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图 11 n=1000, n,.=10 复原 结果 ， 原 图 ( ) 单位 
Jy/pixel, IERI (AE) 单位 Jy/beam, 复原 sd 单位 
Jy/beam, 残 差 (AF) 


Figure 11 For n=1000 and n,=10 deconvolution results: 


Original image (top left, Jy/pixel), dirty image (top right, 
dy/beam), restoration (bottom left, Jy/beam), residual (bot- 
to 
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图 12 n=1000, n,—1. 2. 3. 5 时 复原 结果 (从 上 到 下 ), 复 
原 (Ze) 单位 Jy/beam, 9&25 ( 右 ) 


Figure 12 For n=1000,deconvolution results at n,,—1, 2, 3, and 
5(top to bottom): Restoration (left, Jy/beam), residual (right). 
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SNR of different Reweighted times 
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图 13 图 
迭代 次 数 的 变化 


2.5 
x10? 


(b) 


DR 


2.5 
2.0 
1.5 
1.0 
0.5 


0.0 


13(a) 更 新 权重 次 数 1、2、3、5 次 时 ,SNR PRN ETE, 
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DR of different Reweighted times 
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图 13(b) 更 新 权重 次 数 1、2、3、5 次 时 , DR pë 


Figure 13 In Fig.13(a), the variation of SNR (Signal-to-Noise Ratio) with the number of iterations for 1, 2, 3, and 5 updates of 
weights. In Fig.13(b), the variation of DR (Dynamic Range) with the number of iterations for 1, 2, 3, and 5 updates of weights. 
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14(a) 更 新 权重 次 数 1、2、3、5 次 时 , 在 Dec=—29° 1853” Abit] FTA EE 
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| 线 , 图 14(b) 更 新 权重 次 数 1、2、3、 


5 次 时 , 在 RA=01"20"32° 处 的 剖面 重建 曲线 , 图 14(c) 更 新 权重 次 数 1、2、3、5 次 时 , 在 Dec=—29° 1853” 处 的 章 面 残 差 


曲线 ， 


图 14(d) 更 新 权重 次 数 1、2、3、5 次 时 , 在 RA-0120"32* 处 的 前 面 残 差 曙 线 


Figure 14 Fig.14(a) The section reconstruction curve at Dec=—29°18’53” when the number of update weights is 1, 2, 3, 5 times. 
Fig.14(b) The section reconstruction curve at RA—01/20"32* when the number of update weights is 1, 2, 3, 5 times.Fig.14(c) The 
section residual curve at Dec = —29°18’53” when the number of update weights is 1, 2, 3, 5 times. Fig.14(d) The section residual 
curve at RA—01"20"32* when the number of update weights is 1, 2, 3, 5 times. 
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图 像 的 像素 大 小 为 256*256， 图 像 亮度 值 归 一 
化 到 区 间 [0, 1], 中 心 相 位 为 RA=20°, Dec=-30°, 
观测 使 用 3 个 频率 (100MHz、110MHz 和 120MHz), 
参考 频率 从 100MHz 开始 ， 以 10MHz 的 带宽 递 
增 ， 选 择 频率 通道 1 的 观测 数据 进行 测试 , 结果 
如 图 11 所 示 , 展示 了 M31 的 原始 图 像 、 观 测 12 小 
时 的 脏 图 . 复原 图 像 和 残 差 图 像 . 可 以 看 到 , 天 空 模 
型 与 复原 图 像 基 本 一 致 ， 图 12 展 示 了 不 同 更 新 权重 
次 数 下 的 复原 图 像 和 残留 图 像 对 比 , 在 这 项 测试 中 ， 
迭代 次 数 固定 为 1000 次 ,而 更 新 权重 次 数 分 别 设 
置 为 1、2、3、5 次 . 结果 表明 , 随 着 更 新 权重 次 数 
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越 多 , 对 观测 信号 中 幅度 较 低 的 信号 分 配 的 权重 越 
K, 图 像 的 整体 均 方 根 误差 变 小 , 因此 , 为 了 保障 算 
法 的 收敛 性 和 饱和 性 , 欠 代 次 数 应 该 选 着 合适 的 取 
值 范围 . 为 了 更 直观 地 分 析 不 同 更 新 权重 之 间 的 差 
别 ， 取 复原 图 像 和 残 差 图 像 在 坐标 (120,168)( 真 实 
天 空 图 像 为 RA=01"20"32',Dec=—29° 1853”) 像素 
点 的 剖面 进行 比较 , 如 图 14(a) (b) 所 示 , HE X, 
Y 轴 重 建 图 像 流量 恢复 均 较 好 , 与 原始 天 空 图 像 吻 
合 度 高 . 图 14(c)、(d) 是 剖面 处 X. Y 轴 的 残 差 图 
像 流 量 , 可 以 得 出 更 新 权重 次 数 增加 对 PSF B9573 
抑制 效果 越 好 , 源 的 剩余 流量 已 经 接近 图 像 噪声 水 


的 增加 , 复原 图 像 变 得 更 加 精细 , 残 差 与 噪声 水 平 
趋 于 一 致 . 这 表明 参数 设置 对 于 获得 最 住 射电 图 像 
至 关 重 要 , 并 可 以 根据 具体 需求 进行 调整 以 满足 不 
同 科学 研究 的 要 求 . 在 计算 动态 范围 DR 时 , 通常 
需要 考虑 图 像 峰值 和 计算 选 定 区 域 的 图 像 噪声 , 图 
像 中 的 任何 峰值 像素 值 都 可 视 为 峰值 流量 , 需要 注 
意 的 是 , 图 像 的 分 辩 率 会 影响 到 像素 值 ， 当 分 辨 率 
降低 时 , 每 个 像素 里 包含 流量 增加 ,因为 单个 像素 
覆盖 更 多 的 区 域 , 相对 应 地 , 当 分 辩 率 增加 时 , 相同 
数量 的 流量 分 配 到 了 多 个 像素 上 , 导致 单个 像素 的 
数值 减 小 ， 在 对 256*256 图 像 的 全 图 进行 计算 时 ， 
数值 分 析 结 果 如 图 13(a) 所 示 , 在 固定 更 新 权重 次 
数 为 1、2、3、5 次 的 情况 下 , 迭代 1000 次 可 获得 最 
大 的 信 噪 比 , 达到 23.53dB, 在 迭代 的 初期 , 算法 容 
易 找到 低频 信息 , 从 而 显著 提高 SNR, 然而 , Bee 
代 次 数 的 增加 , 算法 逐渐 双 近 射电 图 像 的 真实 结构 ， 
找到 高 频 信 号 变 得 更 加 困难 , 对 图 像 的 改进 会 变 得 
很 小 , 过 度 增 大 人 迭代 次 数 可 能 会 导致 算法 收敛 性 能 
下 降 , 将 低 于 迭代 阔 值 的 噪声 信号 拟 合 为 射电 源 信 
A, 而 不 是 真实 的 图 像 结 构 . 图 13(b) 显示 , 在 迭代 
1000 次 时 可 获得 最 大 的 动态 范围 , 为 28503.56, 动 
态 范围 受 更 新 权重 次 数 的 影响 较 大 , 更 新 权重 次 数 


x. 


4 总结 与 展望 


为 解决 射电 干涉 成 像 领域 中 的 逆 问 题 , 提出 了 
一 种 基于 压缩 感知 框架 的 射电 图 像 复 原 算 法 ， 与 
传统 的 CLEAN 类 算法 相 比 , 基于 压缩 感知 理论 的 
算法 为 射电 天 文学 领域 的 图 像 重 建 提供 了 一 种 可 借 
鉴 的 方法 , 特别 是 在 处 理 稀 玻 干涉 阵列 观测 数据 时 ， 
利用 目标 射电 源 的 先 验 知识 以 及 在 某 些 域 中 稀 玻 表 
示 来 有 效 克 服 UV 空间 欠 采 样 的 问题 , 从 而 在 射电 
干涉 成 图 中 实现 更 优越 的 重建 效果 . 不 同 的 重 构 算 
法 和 先 验 信息 对 于 图 像 重 建成 功 至 关 重 要 , IUWT- 
FISTA 重建 方法 主要 依赖 于 梯度 下 降 法 , 在 迭代 速 
度 方面 表现 出 色 , 但 易 受 到 噪声 的 干扰 ,其 复原 效 
果 可 能 不 理想 ， 相 比 之 下 ,IUWT-CS 算法 采用 了 
condat 凸 优化 算法 , 能 够 实现 更 加 精细 的 去 噪 和 恢 
复 效 果 . 后 续 研 究 将 着 重 于 开发 适用 于 各 种 射电 源 
的 自 适应 去 卷 积 方法 ,此 外 , 将 研究 如 何在 宽 视 场 
和 大 规模 图 像 的 背景 下 有 效 进行 去 卷 积 处 理 , 以 应 
对 更 复杂 及 更 具 挑战 的 射电 天 文学 场景 . 


SUM 本 研究 使 用 了 中 国 SKA 区 域 中 心 原 型 机 的 资源 进行 大 量 数据 的 模拟 和 测试 . 


参考 文献 


1 Smirnov O M.Revisiting the radio interferometer measurement equation. I. A full-sky Jones formalism.A &A,2011,527:106-118 


xxx-15 


10 


11 


12 


13 


14 


15 


16 


17 


18 


19 


20 


21 


22 


23 


24 


25 


26 


27 


28 


29 
30 


张 讯 等 . 


Smirnov O M.Revisiting the radio interferometer measurement equation. II. Calibration and direction-dependent 
effects. A& A,2011,527:106-118 

Hógbom J A,Cornwell T J. APERTURE SYNTHESIS WITH A NON-REGULAR DISTRIBUTION OF INTERFEROMETER 
BASELINES. Commentary. A&A,1974,500:55-66 

Cornwell T J.Multiscale CLEAN Deconvolution of Radio Synthesis Images.JSTSP,2008,2:793-801 

Haarlem M V.LOFAR: The Low Frequency Array.Eas Publications Series,2005,15:431-444 

Wu X P.White Paper for China SKA (in Chinese). Beijing: Science Press, 2017 [ 武 向 平 . 中 国 SKA 科学 白皮书 . 北京 : 科学 出 
版 社 , 2017] 
An T,Wu X,Lao B,et al.Status and progress of China SKA Regional Centre prototype.Sci Sin-Phys Mech Astron,2022,65 
Wiaux Y,Jacques L,Puy G,et al.Compressed sensing imaging techniques for radio interferometry. MNRAS,2009,395:1733-1742 
Wiaux Y,Puy G, Boursier Y,et al.Spread spectrum for imaging techniques in radio interferometry. MNRAS,2009,400:1029-1038 
Carrillo R E,McEwen J D,Wiaux Y.Sparsity Averaging Reweighted Analysis (SARA): a novel algorithm for ra- 
dio-interferometric imaging. MNRAS,2012,426:1223-1234 

Dabbech  A,Onose A,Abdulaziz  A,et al.Cygnus A  super-resolved via convex optimization from VLA 
data, MNRAS, 2018,476:2853-2866 

Li F,Cornwell T J,Hoog F R.The application of compressive sampling to radio astronomy. I. Deconvolution. A&A,2011,528:10 
Beck A,Teboulle M.A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems.SIIMS,2009,2:183-202 
Condat L.A Primal-Dual Splitting Method for Convex Optimization Involving Lipschitzian, Proximable and Linear Composite 
Terms.J Optim Theory Appl,2013,158:460-479 

Candés E J,Romberg J K,Tao T.Stable signal recovery from incomplete and inaccurate measurements.Commun Pure Appl 
Math,2006,59:1207-1223 


Candés EJ. Compressive sampling. In: Proceedings of the international congress of mathematicians. Madrid,2006,1433-1452 


Candés EJ,Romberg J.Sparsity and incoherence in compressive sampling.Inverse Problems,2006,23:969-985. 

Candés EJ,Wakin M B.An introduction to compressive sampling. IEEE signal processing magazine,2008,25(2):21-30. 

Donoho D L,Huo X.Uncertainty principles and ideal atomic decomposition. [EEE Trans Inf Theory,2001,47:2845-2862 
Thompson A R,Moran J M,Swenson G W. Interferometry and Synthesis in Radio Astronomy.3nd ed.Springer,2017:767-786 
Bhatnagar S,Cornwell T J,Golap K,et al.Correcting direction-dependent gains in the deconvolution of radio interferometric 
images.A& A,2008,487:419-429 

Akiyama K,Ikeda S,Pleau M,et al. Superresolution Full-polarimetric Imaging for Radio Interferometry with Sparse Model- 
ing.The Astronomical Journal,2017:153-159. 

Lu X,Yang J,Yeo T S,et al.Accurate SAR Image Recovery From RFI Contaminated Raw Data by Using Image Domain Mixed 
Regularizations.IEEE Transactions on Geoscience and Remote Sensing,2022,60:1-13. 

Starck J,Murtagh F,Bertero M.Starlet transform in astronomical data processing.New York:Springer,2015:2053-2098 

Candés E J,Recht B.Exact Matrix Completion via Convex Optimization.FoCM,2009,9:717-772 

McMullin J P,Waters B R,Schiebel D,et al, CASA Architecture and Applications.In: Astronomical Data Analysis Software and 
Systems XVI ASP Conference Series. Arizona,2007 

Guo $ G, Lu Y, An T, et al.Scientific data flow and array simulation analysis for the SKA1 era (in Chinese).Sci Sin-Phys 
Mech Astron,2023, 53: 229504,ChinaXiv: 202206.00184 [3824 56, 陆 扬 , 安 涛 , 等 . 面向 SKA1 时 代 的 科学 数据 流 及 阵列 模拟 分 析 . 
国 科学 : 物理 学 力学 天 文学 ,2023, 53: 229504, ChinaXiv: 202206.00184] 

Li F,Brown S,Cornwell T J,et al.The application of compressive sampling to radio astronomy. II. Faraday rotation measure 
synthesis. A&A ,2011,531:8 

Briggs D S.High fidelity deconvolution of moderately resolved sources. The New Mexico Institute of Mining and Technology,1995 
Dabbech A,Ferrari C,Mary D L,et al. MORESANE: MOdel REconstruction by Synthesis-ANalysis Estimators-A sparse decon- 
volution algorithm for radio interferometric imaging. A&A,2014,576:16 


g 


O 


xxx-16 


张 讯 等 . 


A radio astronomy image restoration algorithm based on 
compressed sensing framework 
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Deconvolution of radio astronomy images is a key data processing technique. Its main goal is to remove 
the effects introduced by the instrument from the observed sky images to recover the original sky images. 
However, radio interferometer arrays employ sparse interferometric arrays, whose imaging principles differ 
from those of optical telescopes. If the sampling points in the UV space are not sufficiently dense, this will 
lead to insufficient high-resolution information during image reconstruction. Traditional radio astronomy 
image reconstruction algorithms fail to fundamentally solve the problem of UV space undersampling. This 
paper adopts the compressed sensing theoretical framework, combines prior knowledge of the sparsity of 
radio astronomy images, and studies a new radio astronomy image deconvolution algorithm, namely the 
IUW'T-CS-based radio interferometric image reconstruction algorithm. This algorithm transforms the dirty 
image deconvolution process into a convex optimization problem to find the global minimum. To evaluate 
the reconstruction performance of this algorithm, we used the OSKAR radio astronomy simulation software 
package to simulate low SKA1 observations and performed deconvolution processing on the extended radio 
sources obtained. Experimental results show that, compared with the HOGBOM-CLEAN, MS-CLEAN, 
and IUWT-FISTA methods, the IUWT-CS method significantly improves the reconstruction quality of radio 
images and achieves finer denoising and restoration effects. 


image deconvolution, compressed sensing, sparse representation, image reconstruction 
PACS: 95.75.Mn,07.05.Pj,95.75.Tv 
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